2 Introduction

In the workshop, we learned about:

• What is Spatial Data?
• What is the sf framework for R?

To delve in further, let’s see some spatial data in action.

We’ll work with the sf library first.

library(sf)

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

2.1 Load Spatial Data

Chi_tracts = st_read("data/geo_export_aae47441-adab-4aca-8cb0-2e0c0114096e.shp")

## Reading layer `geo_export_aae47441-adab-4aca-8cb0-2e0c0114096e' from data source 
##   `/Users/maryniakolak/code/Intro2RSpatialMed/data/geo_export_aae47441-adab-4aca-8cb0-2e0c0114096e.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 801 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.94025 ymin: 41.64429 xmax: -87.52366 ymax: 42.02392
## CRS:           4326

2.2 Non-Spatial & Spatial Views

head(Chi_tracts)

## Simple feature collection with 6 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.68822 ymin: 41.72902 xmax: -87.62394 ymax: 41.87455
## CRS:           4326
##   commarea commarea_n countyfp10     geoid10 name10        namelsad10 notes
## 1       44         44        031 17031842400   8424 Census Tract 8424  <NA>
## 2       59         59        031 17031840300   8403 Census Tract 8403  <NA>
## 3       34         34        031 17031841100   8411 Census Tract 8411  <NA>
## 4       31         31        031 17031841200   8412 Census Tract 8412  <NA>
## 5       32         32        031 17031839000   8390 Census Tract 8390  <NA>
## 6       28         28        031 17031838200   8382 Census Tract 8382  <NA>
##   statefp10 tractce10                       geometry
## 1        17    842400 POLYGON ((-87.62405 41.7302...
## 2        17    840300 POLYGON ((-87.68608 41.8229...
## 3        17    841100 POLYGON ((-87.62935 41.8528...
## 4        17    841200 POLYGON ((-87.68813 41.8556...
## 5        17    839000 POLYGON ((-87.63312 41.8744...
## 6        17    838200 POLYGON ((-87.66782 41.8741...

plot(Chi_tracts)

2.3 Spatial Data Structure

str(Chi_tracts)

## Classes 'sf' and 'data.frame':   801 obs. of  10 variables:
##  $ commarea  : Factor w/ 77 levels "1","10","11",..: 39 55 28 25 26 21 62 49 74 75 ...
##  $ commarea_n: num  44 59 34 31 32 28 65 53 76 77 ...
##  $ countyfp10: Factor w/ 1 level "031": 1 1 1 1 1 1 1 1 1 1 ...
##  $ geoid10   : Factor w/ 801 levels "17031010100",..: 785 767 772 773 756 751 584 513 684 34 ...
##  $ name10    : Factor w/ 801 levels "1001","1002",..: 782 764 769 770 753 748 545 443 663 266 ...
##  $ namelsad10: Factor w/ 801 levels "Census Tract 1001",..: 782 764 769 770 753 748 545 443 663 266 ...
##  $ notes     : Factor w/ 7 levels "Half in CA 64 (Midway Airport)",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ statefp10 : Factor w/ 1 level "17": 1 1 1 1 1 1 1 1 1 1 ...
##  $ tractce10 : Factor w/ 801 levels "010100","010201",..: 785 767 772 773 756 751 584 513 684 34 ...
##  $ geometry  :sfc_POLYGON of length 801; first list element: List of 1
##   ..$ : num [1:243, 1:2] -87.6 -87.6 -87.6 -87.6 -87.6 ...
##   ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA
##   ..- attr(*, "names")= chr  "commarea" "commarea_n" "countyfp10" "geoid10" ...

st_crs(Chi_tracts)

## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS84(DD)",
##     DATUM["WGS84",
##         SPHEROID["WGS84",6378137.0,298.257223563]],
##     PRIMEM["Greenwich",0.0],
##     UNIT["degree",0.017453292519943295],
##     AXIS["Geodetic longitude",EAST],
##     AXIS["Geodetic latitude",NORTH]]

2.4 Exploring Coordinate Reference Systems

Chi_tracts.moll <- st_transform(Chi_tracts, crs = "+proj=moll")
plot(st_geometry(Chi_tracts.moll), border = "gray", lwd = 2, main = "Mollweide", sub="preserves areas")

st_crs(Chi_tracts.moll)

## Coordinate Reference System:
##   User input: +proj=moll 
##   wkt:
## PROJCS["unnamed",
##     GEOGCS["WGS 84",
##         DATUM["unknown",
##             SPHEROID["WGS84",6378137,298.257223563]],
##         PRIMEM["Greenwich",0],
##         UNIT["degree",0.0174532925199433]],
##     PROJECTION["Mollweide"],
##     PARAMETER["central_meridian",0],
##     PARAMETER["false_easting",0],
##     PARAMETER["false_northing",0]]

Chi_tracts.54019 = st_transform(Chi_tracts, 54019)
plot(st_geometry(Chi_tracts.54019), border = "gray", lwd = 2, main = "Winkel", sub="minimal distortion")

st_crs(Chi_tracts.54019)

## Coordinate Reference System:
##   User input: EPSG:54019 
##   wkt:
## PROJCS["World_Winkel_II",
##     GEOGCS["GCS_WGS_1984",
##         DATUM["WGS_1984",
##             SPHEROID["WGS_84",6378137.0,298.257223563]],
##         PRIMEM["Greenwich",0.0],
##         UNIT["Degree",0.0174532925199433]],
##     PROJECTION["Winkel_II"],
##     PARAMETER["False_Easting",0.0],
##     PARAMETER["False_Northing",0.0],
##     PARAMETER["Central_Meridian",0.0],
##     PARAMETER["Standard_Parallel_1",50.45977625218981],
##     UNIT["Meter",1.0],
##     AUTHORITY["Esri","54019"]]

Chi_tracts.3435 <- st_transform(Chi_tracts, 3435)
st_crs(Chi_tracts.3435)

## Coordinate Reference System:
##   User input: EPSG:3435 
##   wkt:
## PROJCS["NAD83 / Illinois East (ftUS)",
##     GEOGCS["NAD83",
##         DATUM["North_American_Datum_1983",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6269"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4269"]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",36.66666666666666],
##     PARAMETER["central_meridian",-88.33333333333333],
##     PARAMETER["scale_factor",0.999975],
##     PARAMETER["false_easting",984250.0000000002],
##     PARAMETER["false_northing",0],
##     UNIT["US survey foot",0.3048006096012192,
##         AUTHORITY["EPSG","9003"]],
##     AXIS["X",EAST],
##     AXIS["Y",NORTH],
##     AUTHORITY["EPSG","3435"]]

plot(st_geometry(Chi_tracts.3435), border = "gray", lwd = 2, main = "NAD83 / Illinois East (ftUS)", sub="topo mapping & survey use")

Chi_tracts.Hawaii = st_transform(Chi_tracts, 102114)
plot(st_geometry(Chi_tracts.Hawaii), border = "gray", lwd = 2, main = "Old Hawaiian UTM Zone 4N", sub="wrong projection!")

2.5 Refine Basic Map

library(tmap)

## Warning: replacing previous import 'sf::st_make_valid' by
## 'lwgeom::st_make_valid' when loading 'tmap'

tm_shape(Chi_tracts) + tm_borders(alpha=0.5) 

tm_shape(Chi_tracts) + tm_fill(col = "gray90") + tm_borders(alpha=0.2, col = "gray10") +
  tm_scale_bar(position = ("left"), lwd = 0.8) +
  tm_layout(frame = F)

Check out https://rdrr.io/cran/tmap/man/tm_polygons.html for more ideas.

2.6 Arrange multiple maps

tracts.4326 <- tm_shape(Chi_tracts) + tm_fill(col = "gray90") +
  tm_layout(frame = F, title = "EPSG 4326")

tracts.54019 <- tm_shape(Chi_tracts.54019) + tm_fill(col = "gray90") +  tm_layout(frame = F, title = "EPSG 54019")

tmap_arrange(tracts.4326, tracts.54019)

2.7 Interactive Mode

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(Chi_tracts) + tm_fill(col = "gray90") + tm_borders(alpha=0.2, col = "gray10") +
  tm_scale_bar(position = ("left"), lwd = 0.8) +
  tm_layout(frame = F)

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string,
## which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is
## deprecated.

tm_shape(Chi_tracts) + tm_fill(col = "gray90", alpha = 0.5) + tm_borders(alpha=0.2, col = "gray10") + tm_scale_bar(position = ("left"), lwd = 0.8) +  tm_layout(frame = F)

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string,
## which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is
## deprecated.

tmap_mode("plot")

## tmap mode set to plotting

2.8 Overlay Zip Code Boundaries

First, we read in zip code boundaries. This data was downloaded directly from the City of Chicago Data Portal as a shapefile.

Chi_Zips = st_read("data/geo_export_54bc15d8-5ef5-40e4-8f72-bb0c6dbac9a5.shp")

## Reading layer `geo_export_54bc15d8-5ef5-40e4-8f72-bb0c6dbac9a5' from data source 
##   `/Users/maryniakolak/code/Intro2RSpatialMed/data/geo_export_54bc15d8-5ef5-40e4-8f72-bb0c6dbac9a5.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 61 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
## CRS:           4326

Next, we layer the new shape in – on top of the tracts. We use a thicker border, and try out a new color. Experiment!

## FIRST LAYER: CENSUS TRACT BOUNADRIES
tm_shape(Chi_tracts.3435) + tm_fill(col = "gray90") +
  tm_borders(alpha=0.2, col = "gray10") + 

## SECOND LAYER: ZIP CODE BOUNDARIES WITH LABEL
tm_shape(Chi_Zips) + tm_borders(lwd = 2, col = "#0099CC") +
  tm_text("zip", size = 0.7) +
  
## MORE CARTOGRAPHIC STYLE
  tm_scale_bar(position = ("left"), lwd = 0.8) +
  tm_layout(frame = F)

More Resources

On spatial data basics & sf:

• https://geocompr.robinlovelace.net/intro.html

• https://geodacenter.github.io/opioid-environment-toolkit/spatial-data-introduction.html

On projections:

• https://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/projection-basics-for-gis-professionals.htm

• https://geocompr.robinlovelace.net/reproj-geo-data.html

• https://datacarpentry.org/organization-geospatial/03-crs/index.html

On tmap:

• https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html

• https://geocompr.robinlovelace.net/adv-map.html